Vortex glass transition in a frustrated 3D XY model with disorder 
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The anisotropic frustrated three dimensional (3D) XY model with disorder in the coupling con- 
stants is simulated as a model of a point disordered superconductor in an applied magnetic field. 
From a finite size scaling analysis of the helicity modulus it is concluded that the data is consistent 
with a finite temperature transition with isotropic scaling and the correlation length exponent is 
found to be v — 1 .50 ± 0.f2, consistent with 3D gauge glass universality. 
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The hypothesis of a vortex glass in disordered high 
temperature superconductors \& Q has spurred much 
research and many discussions during more than one 
decade and continues to be a very controversial issueQ. 
The essence of this suggestion is that random point disor- 
der in superconductors may conspire with the vortex line 
interaction to pin the vortices and that this takes place 
through a sharp transition into a phase with vanishing 
linear resistance. 

Computer simulations have for quite some time played 
an important role in the examination of critical phenom- 
ena and have also recently been used in the study of some 
vortex glass models. One important such model is the 
three dimensional (3D) gauge glass, which is an isotropic 
3D XY model with randomness included through a ran- 
dom vector potential added to the phase difference of the 
superconducting order parameter. The evidence has for 
quite some time pointed at a finite temperature transition 
in this model but strong evidence for a real phase 

transition has been obtained only recently^ through the 
use of the exchange Monte Carlo (MC) technique|7]. The 
value of the correlation length exponent was then found 
to be v — 1.39 ± 0.20. It is however generally recognized 
that the 3D gauge glass models is too much of a simplifi- 
cation to allow for any safe conclusions regarding the be- 
havior of disordered superconductors in applied magnetic 
fields Most seriously, the gauge glass is an isotropic 
model with no net field which means that the possibility 
of anisotropic scaling is excluded at the outset. 

Two studies with the necessary ingredients of disorder 
and applied field have so far been reported in the litera- 
ture. The first is an examination of a frustrated 3D XY 
model with randomness in the couplings 0, and the data 
was there interpreted as evidence for a phase transition 
with v w 2.2. The crossing of the data for different sizes 
expected from finite size scaling, was however not en- 
tirely convincing, possibly because of the open boundary 
conditions employed in the simulation. The second study 
gives results for a 3D random pinning model with strong 
disorder 9] . The value of the correlation length exponent 
was there found to be v ss 0.7, undistinguishable from 
v w 0.67 in the pure zero field 3D XY model. This is at 
odds with the common expectation that a vortex glass 



transition should be in a different universality class than 
the pure model. 

In this paper we present results from large scale sim- 
ulations on a frustrated 3D XY model with disorder 
in the coupling constants. The quantity in focus is 
the helicity modulus and we find that the data is con- 
sistent with a finite temperature glass transition with 
isotropic scaling and obtain the correlation length expo- 
nent v = 1.50±0.12. The agreement with v = 1.39±0.20 
for the 3D gauge glass model 6] suggests a common uni- 
versality class. 

The model we simulate is given by the Hamiltonian|l0| 

H = - ^2 J v cos (^ _ _ A ^ + <W> U) 

bonds i^i 

where 9i is the phase of the superconducting wave func- 
tion at site i of a periodic L x x L y x L z lattice and the 
sum is over all bonds in directions fj, = x, y, z. An ap- 
plied magnetic field in the z direction corresponding to 
1/5 flux quantum per plaquette is obtained through the 
quenched vector potential with the choice Ai X = yi 2tt/5, 
and Ai y = Ai Z = 0. The randomness is included through 
disorder in the coupling constants, 

Jin = J±(l + pe^), /i = x, y, 

where ti^ are independent variables from a Gaussian 
distribution with (e^) = and (e^) 2 = 1. The dis- 
order strength p — 0.4 together with the anisotropy 
Jii = Ji/40 were choosen since they were found suf- 
ficient to prohibit the formation of Abrikosov lattices. 
The simulations are performed with fluctuating twist 
boundary conditions [111 ; the 5^ in Eq. ifjj are the twist 
variables and the total twist in the respective directions 
are = L^5^. The simulations are performed with 
L = L x = L y and a fixed aspect ratio, L/ L z — 5/3. The 
temperature is given in units of Jj_. 

The quantity in focus in our analysis is the helicity 
modulus which is defined from the change in free energy 
density, /, due to an applied twist, <5 M : T p = d 2 f/dS 2 
[l2T |. To use the helicity modulus as a signal of the stiff- 
ness of the system the derivative should be evaluated 
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at the twist that minimizes the free energy. In ordered 
systems this minimum is always at zero twist and the 
helicity modulus may then be evaluated by means of a 
correlation function determined with periodic boundary 
conditions, A M = 0. For disordered systems, however, 
the minimizing twist will in general be different from zero 
and one then needs to make simulations with the twist 
variables A M as additional dynamical variables and col- 
lect histograms P A1 (A A1 ). The helicity modulii are then 
determined from the free energies P M = — TlnP M , as dis- 
cussed below. To analyze the critical behavior we use the 
standard scaling relation for the helicity modulus in 3D, 

LT ~ g (tL 1 '") , (2) 

where t — (T— T c ) /T c and g is a scaling functional. This 
expression presumes isotropic scaling; the more general 
scaling relations are gives in Ref. j^J- A naive analysis 
of T(L) rather than the correct scaling quantity LT, led 
to an erroneous conclusion regarding the existence of a 
vortex glass phase in Ref. 01 • 
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TABLE I: Parameters describing the simulations. For systems 
of size LxLxLz we simulated Nd disorder configurations with 
Nt temperatures in the range T m i n < T < T max , cf. Eq. (|3J. 
The acceptance ratio for the exchange step is given by X acc . 
Of the bins corresponding to 2 18 = 262144 sweeps r cq are 
first discarded and the remaining r ma x are used for calculating 
averages. The same information is also given in terms of the 
number of sweeps for equilibration and for collecting data. 

Our simulations are performed with exchange MC 
which is a method for simultaneously simulating mul- 
tiple copies of a particular configuration of disorder with 
each copy at a different temperature. According to cer- 
tain rules these copies may now and then interchange 
temperature^ and therefore effectively perform random 
walks in temperature space. These random changes in 
temperature greatly help the different copies avoid get- 
ting trapped in restricted parts of the phase space and 
therefore makes it possible to sample the whole phase 
space and obtain the true thermodynamic averages. In 
our simulations the temperatures were chosen according 
to the equation 

(rp \ m/Nr 
, m = 0,...,iV T -l, (3) 
min / 

with T ma x = 0.24 and T m i n as given in Table [I] Before 
doing the actual exchange MC the initial spin configura- 
tions for the Nt temperatures were obtained by slowly 



cooling the system with standard MC simulations. The 
number of temperatures, the acceptance ratio for the ex- 
change step, and the number of disorders simulated for 
the different system sizes are shown in Table [|] The ex- 
change steps are attempted once every 16 sweep. 

The exchange MC method ensures that all the differ- 
ent copies remain at thermal equilibrium as soon as equi- 
librium has been reached. The approach to equilibrium 
may however be very slow since information and con- 
figurations have to propagate all the way from high to 
low temperatures. We have carefully examined the ap- 
proach to equilibrium and especially for our largest sizes 
the times for equilibration are indeed very long. For the 
next largest size, L = 20, equilibration is only reached 
after about 2.9 x 10 6 sweeps, and to make thermalization 
at all possible for our largest size, L = 25, we chose not 
to go to quite that low temperatures for the largest size, 
cf. Table [Q This was decided since the time required for 
thermalization may increase very rapidly with decreasing 
temperature. Still, about 4.5 x 10 6 sweeps were necessary 
to reach equilibrium for L = 25. 
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FIG. 1: P X (A X ) for a certain disorder realisation. Note that 
the peaks of the histograms become higher and sharper as 
the temperature is lowered. This kind of data is used for the 
determination of T^ through Eq. 10. The inset gives the 
corresponding quantity for the z direction. 

For each system size, disorder configuration and tem- 
perature the main output from the simulations are his- 
tograms P M (A M ; t), where r enumerates bins correspond- 
ing to 2 18 = 262144 sweeps over the lattice. The further 
analysis is then based on the average 

iyA M H-^ T £ X p M (A M ;T), (4) 

'max i 

T— 1 

and the helicity modulus is determined from the curva- 
ture at the minimum of the associated free energy by 
fitting a second order polynomial to the free energy in a 
narrow interval around the minimum, A°, which is also 
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determined in the fit, 
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However, it turns out that the values of T M obtained in 
this way are biased towards too large values. The origin 
of this bias as well as the method employed to eliminate 
it from the data is discussed below after the discussion 
of the results. 

We now focus on disorder averaged quantities for which 
the bias mentioned above has already been eliminated. 
With [. . .] av denoting the average over disorder configu- 
rations we define 
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Our results for LT^ are shown in Fig. EI T° a very good 
accuracy the data for the different sizes cross at a single 
temperature. To further verify the scaling according to 
Eq. (J2J we fit our data for LT± near T c to a fourth order 
polynomial expansion of g(tL x ' v ) and obtain the values 
v = 1.50 ± 0.12 and T c = 0.123 ± 0.003. The collapse 
which is shown in Fig. [3] is excellent and holds in a sur- 
prisingly large temperature interval. The error estimates 
are obtained with a resampling technique and correspond 
to one standard deviation. 
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FIG. 2: LT± versus temperature for four different system 
sizes. The curves cross at a single point, which is an indication 
of critical behavior. The inset shows LTn/Jn which shows a 
similar crossing at almost the same temperature. 



A vortex glass transition should also be seen in the 
parallel component of the helicity modulus and since the 
anisotropy exponent enters the scaling relation in differ- 
ent ways for T|| and T^ 9|, a scaling analysis of £Tii 
constitutes an additional test that the scaling actually is 
isotropic. The crossing of LTii at T c for different sizes 
and the scaling collapse are shown in the insets of Figs. El 
and |3 respectively. From the not so smooth curves it is 
clear that the precision in T|| is much worse than for Tj_. 
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FIG. 3: The scaling collapse of LT± gives v = 1.50 ± 0.12 
and T c = 0.123 ± 0.003. The inset is an attempt to collapse 
LT||/J|| with the same values for v and T c . The nice col- 
lapse gives additional support for a vortex glass transition 
with isotropic scaling. 



This may be traced back to the less good quality of the 
histograms P Z (A Z ), as shown in the inset of Fig.Q] Nev- 
ertheless, the data for LTii (see inset in Fig- collapses 
nicely when using the same values of v and T c as in the 
main figure. This therefore constitutes additional evi- 
dence for a vortex glass transition in the model. We now 
turn to the bias mentioned above and also discuss simi- 
larities and differences with the more common method to 
determine the root mean square current / rms , in analyses 
of models with disorder. 

As mentioned above the determination of T from a 
twist histogram suffers from a bias towards too large val- 
ues. To investigate the reason for this bias we performed 
additional simulations with fluctuating twist boundary 
conditions of the pure zero-field 3D XY model and col- 
lected histgrams -P(A, r) where r enumerates the bins. 
This data was then used to calculate averages over r avcr 
consecutive bins -P(A; r avcr ), which in turn were used to 
determine the helicity modulus. By repeating this pro- 
cedure for several values of T avcr the dependence of T on 
Tavcr was determined. From this kind of analysis it was 
found that the bias decays as 1/ T ave r to a very good pre- 
cision. Another finding is that the bias may be made to 
vanish altogether by making use of A = in Eq. (J5J 
instead of using A as a free parameter. (A = is 
the known value of the minimizing twist in the pure sys- 
tem.) The latter observation suggests that the bias is re- 
lated to the usual complication in determining the width 
(variance) of a distribution when the true average is not 
known 14J. Since T is inversely related to the variance 
of the distribution -P(A) around the maximum it follows 
that the estimates based on runs of length r a y Br would be 
expected to decay towards the true value as|l4j 
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where 6 is a free parameter related to the decorrelation 
time in the simulations. For small values of 6/r avcr the 
expected behavior is therefore entirely in accordance with 
the l/r aver decay discussed above. 

Returning to our data for the vortex glass model, the 
procedure used to determine the data points in Fig. |3 
consists of three steps: (i) Determine T M (r avcr ) for each 
disorder configuration and several values of r avor by fit- 
ting histogram P M (A M ; r avC r) based on T avcr consecutive 
bins, P m (A m ,t) to Eq. (0). (ii) Calculate the disorder 
averaged quantities T^(T ave r) and Ti|(r aver ), cf. Eqs. 10. 
(iii) Fit this data to Eq. (J7J) to obtain the unbiased es- 
timates Tj_ = T^(oo) and Ty = T||(oo). The last step 
is illustrated in Fig. |1 for T = 0.125 close to T c . The 
error bars on the last point for each size are the errors 
associated with the disorder average. 
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Note also that there is a bias in the determination of 
the root mean square current / rms that is very similar to 
the bias in T discussed above. For I rms = ^/[J 2 ] av the 
origin of this bias is that the statistical error SI gives a 
term (SI) 2 > 0. Since SI would be expected to vanish 
with simulation time as ~ 1 /^/r avcr the bias in J rms would 
vanish as l/r aver , which is essentially the same as the 
behavior of T(T avC r). There is however a difference in 
that the bias in J rms is easily eliminated by performing 
two independent simulations, a and 3, for each disorder 
and measuring the quantity [/ Q / ) 3] a v|6|. Since the bias 
in T is of a very different origin it cannot be eliminated 
with such methods. 

To summarize, we have performed a finite size scaling 
analysis of the helicity modulus in a frustrated 3D XY 
model with disorder in the coupling constants. The data 
is consistent with isotropic scaling and the correlation 
length exponent is found to be v = 1.50 ± 0.12. The 
good agreement with v = 1.39 ± 0.20 of the 3D gauge 
glass suggests that the two models actually do belong to 
same universality class. 
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FIG. 4: The figure shows the elimination of the bias in Ti 
by fitting Ti(T avcr ) to Eq. with b and Tx(oo) as free 
parameters. The time is in units of the bin size which is 
2 18 = 262144 sweeps. 

We now shortly discuss the relation between the helic- 
ity modulus calculated in the present paper and the more 
common method to determine the root mean square cur- 
rent / rms [j|. Since the current and the helicity modu- 
lus are first and second derivatives, respectively, of the 
same function F(A), both quantities effectively probe the 
roughness of the function F(A). The low temperature 
phase is characterized by large energy barriers growing 
with system size whereas the free energy above the tran- 
sition temperature becomes a flat function of A in the 
limit of large L. However, T turns out to be more effi- 
cient in measuring the roughness of F(A). The reason is 
that T measures a property at an extremum (the mini- 
mum of the free energy) whereas the current at A = 
for a given shape of the function may be large or small 
depending on the location of the structure in F(A). I rms 
is therefore only on the average a good measure of the 
properties of F(A) and this has to be compensated for 
by using a larger number of disorder configurations. This 
is the reason for the good precision in our data in spite 
of the rather small number of disorder configurations. 
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